package scu.maqiang.numeric;

import java.util.Arrays;

import scu.maqiang.mesh.Tecplot;

public class FourierTransformation {

	public static void main(String[] args) {
/*		double[] x1 = {1, 0, -1, 0};
		Complex[] result1 = DFT(x1);
		System.out.println(Complex.toString(result1));
		double[] x2 = {0, 1, 0, -1};
		Complex[] result2 = DFT(x2);
		System.out.println(Complex.toString(result2));
		double[] x3 = {1, 1, 1, 1};
		Complex[] result3 = DFT(x3);
		System.out.println(Complex.toString(result3));
		double[] x4 = {0,-1, 0, 1};
		Complex[] result4 = DFT(x4);
		System.out.println(Complex.toString(result4));
		double[] x5 = {0, 1, 0, -1, 0, 1, 0, -1};
		Complex[] result5 = DFT(x5);
		System.out.println(Complex.toString(result5));
		
		Complex[] r1 = new Complex[4];
		r1[0] = new Complex(1, 0);
		r1[1] = new Complex(0, 0);
		r1[2] = new Complex(0, 0);
		r1[3] = new Complex(0, 0);
		Complex[] rr1 = inverseDFT(r1);
		System.out.println(Complex.toString(rr1));
		
		r1[1].real = 1.0;
		r1[2].real = -1.0;
		r1[3].real = 1.0;
		Complex[] rr2 = inverseDFT(r1);
		System.out.println(Complex.toString(rr2));
		r1[1] = new Complex(0.0, -1);
		r1[2].real = 1.0;
		r1[3] = new Complex(0.0, 1.0);
		Complex[] rr3 = inverseDFT(r1);
		System.out.println(Complex.toString(rr3));
		
		Complex[] r2 = new Complex[8];
		r2[0] = new Complex(1.0, 0.0);
		r2[1] = new Complex(0.0, 0.0);
		r2[2] = new Complex(-0.5, 0.0);
		r2[3] = new Complex(0.0, 0.0);
		r2[4] = new Complex(1.0, 0.0);
		r2[5] = new Complex(0.0, 0.0);
		r2[6] = new Complex(-0.5, 0.0);
		r2[7] = new Complex(0.0, 0.0);
		Complex[] rr4 = inverseDFT(r2);
		System.out.println(Complex.toString(rr4));
		
		rr1 = new Complex[4];
		rr1[0] = new Complex(1, 0);
		rr1[1] = new Complex(0, 0);
		rr1[2] = new Complex(-1, 0);
		rr1[3] = new Complex(0, 0);
		result1 = DFT(rr1);
		System.out.println(Complex.toString(result1));
		
		double[] case103 = {-2.2, -2.8, -6.1, -3.9, 0.0, 1.1, -0.6, -1.1};
		result1 = DFT(case103);
		System.out.println(Complex.toString(result1));
		
		double[] t = MVO.linspace(0, 1, 129);
		double[] interY = new double[129];
		double[] inter = {0, 1};
		Complex[] myResult = result1;
		Arrays.setAll(interY, i -> TrigonometricInterpolation(myResult, inter, t[i]));
		Tecplot.LineXY("TrigonInterResult.dat", t, new double[][] {interY});*/

		/*
		double[][] xy = {{1, 0}, {0, 0}, {-1, 0}, {0, 0}};
		Complex[] source = Complex.fromTwoArrays(xy);
		Complex[] result = fft(source);
		System.out.println(Complex.toString(result));
		
		double[][] case103 = {{-2.2,0}, {-2.8, 0.0}, {-6.1, 0.0}, {-3.9, 0.0}, {0.0, 0.0}, 
				              {1.1, 0.0}, {-0.6, 0.0}, {-1.1, 0.0}};
		source = Complex.fromTwoArrays(case103);
		result = fft(source);
		System.out.println(Complex.toString(result));
		source = ifft(result);
		System.out.println(Complex.toString(source));*/
		
		double[][] x = {{0.709364830858073, 0.655098003973841, 0.959743958516081, 0.751267059305653}, 
				   {0.754686681982361, 0.162611735194631, 0.340385726666133, 0.255095115459269}, 
				   {0.276025076998578, 0.118997681558377, 0.585267750979777, 0.505957051665142}, 
				   {0.679702676853675, 0.498364051982143, 0.223811939491137, 0.699076722656686}};
		Complex[][] xx = new Complex[4][4];
		for(int i = 0; i < 4; i++) {
			for(int j = 0; j < 4; j++) {
				xx[i][j] = new Complex(x[i][j], 0.0);
			}
		}
		
		Complex[][] ffxx = ifft2(xx);
		for(int i = 0; i < 4; i++) {
			for(int j = 0; j < 4; j++) {
				System.out.print(ffxx[i][j] + "\t");
			}
			System.out.println();
		}
		
	}
	
	public static Complex[] DFT(double[] x) {
		int n = x.length; //设n是偶数
		Complex[] result = new Complex[n];
		for(int i = 0; i < n; i++) {
			result[i] = new Complex(0, 0);
		}
		result[0].real = MVO.sum(x)/ Math.sqrt(n);
		result[0].image = 0.0;
		for(int i = 1; i <= n/2; i++) {
			result[i].real = 0.0;
			result[i].image = 0.0;
			for(int j = 0; j < n; j++) {
				result[i].real += x[j] * Math.cos(2 * j * i * Math.PI / n);
				result[i].image -= x[j] * Math.sin(2 * j * i * Math.PI / n);
			}
			result[i].real /= Math.sqrt(n);
			result[i].image /= Math.sqrt(n);
			result[n - i].real = result[i].real;
			result[n - i].image = -result[i].image;
		}
		return result;
	}
	
	public static Complex[] DFT(Complex[] x) {
		int n = x.length; //设n是偶数
		Complex[] result = new Complex[n];
		for(int i = 0; i < n; i++) {
			result[i] = new Complex(0, 0);
		}
		Pair<double[], double[]> xir = Complex.toTwoArrays(x);
		double[] xi = xir.getFirst();
		double[] xr = xir.getSecond();
		result[0].real = MVO.sum(xi) / Math.sqrt(n);
		result[0].image = MVO.sum(xr) / Math.sqrt(n);
		double wr, wi;
		for(int i = 1; i < n; i++) {
			result[i].real = 0.0;
			result[i].image = 0.0;
			for(int j = 0; j < n; j++) {
				wr = Math.cos(2 * j * i * Math.PI / n);
				wi = -Math.sin(2 * j * i * Math.PI / n);
				result[i].real += xi[j] * wr - xr[j] * wi;
				result[i].image += xi[j] * wi + xr[j] * wr;
			}
			result[i].real /= Math.sqrt(n);
			result[i].image /= Math.sqrt(n);
		}
		return result;
	}
	
	public static Complex[] inverseDFT(Complex[] x) {
		int n = x.length; //设n是偶数
		Complex[] result = new Complex[n];
		for(int i = 0; i < n; i++) {
			result[i] = new Complex(0, 0);
		}
		Pair<double[], double[]> xir = Complex.toTwoArrays(x);
		double[] xi = xir.getFirst();
		double[] xr = xir.getSecond();
		result[0].real = MVO.sum(xi) / Math.sqrt(n);
		result[0].image = MVO.sum(xr) / Math.sqrt(n);
		double wr, wi;
		for(int i = 1; i < n; i++) {
			result[i].real = 0.0;
			result[i].image = 0.0;
			for(int j = 0; j < n; j++) {
				wr = Math.cos(2 * j * i * Math.PI / n);
				wi = Math.sin(2 * j * i * Math.PI / n);
				result[i].real += xi[j] * wr - xr[j] * wi;
				result[i].image += xi[j] * wi + xr[j] * wr;
			}
			result[i].real /= Math.sqrt(n);
			result[i].image /= Math.sqrt(n);
		}
		return result;
	}
	
	public static double TrigonometricInterpolation(Complex[] y, double[] interval, double t) {
		double c = interval[0];
		double d = interval[1];
		double position = 2 * Math.PI * (t - c) / (d - c);
		double result = 0.0;
		int n = y.length;
		for(int i = 0; i < n; i++) {
			result += y[i].real * Math.cos(i * position);
			result -= y[i].image * Math.sin(i * position);
		}
		return result / Math.sqrt(n);
	}
	
	/**
	 * 使用快速Fourier变换算法对复数向量作离散Fourier变换
	 * @param source 复数向量, 要求长度n必须为2的幂次, 即n = 2^k, k为整数
	 * @return 离散Fourier变换后的向量
	 */
	public static Complex[] fft(Complex[] source) {
		int i, j, k, n1, n2, a;
		double c, s, t1, t2;
		int n = source.length;
		Pair<double[], double[]> twoArrays = Complex.toTwoArrays(source);
		double[] x = twoArrays.getFirst();
		double[] y = twoArrays.getSecond();
		
		int m = (int) (Math.log(n) / Math.log(2));

		// Make sure n is a power of 2
		if (n != (1 << m))
			throw new RuntimeException("FFT length must be power of 2");

		// precompute tables
		double[] cos = new double[n / 2];
		double[] sin = new double[n / 2];

		for (i = 0; i < n / 2; i++) {
			cos[i] = Math.cos(-2 * Math.PI * i / n);
			sin[i] = Math.sin(-2 * Math.PI * i / n);
		}
		// Bit-reverse
		j = 0;
		n2 = n / 2;
		for (i = 1; i < n - 1; i++) {
			n1 = n2;
			while (j >= n1) {
				j = j - n1;
				n1 = n1 / 2;
			}
			j = j + n1;

			if (i < j) {
				t1 = x[i];
				x[i] = x[j];
				x[j] = t1;
				t1 = y[i];
				y[i] = y[j];
				y[j] = t1;
			}
		}

		// FFT
		n1 = 0;
		n2 = 1;

		for (i = 0; i < m; i++) {
			n1 = n2;
			n2 = n2 + n2;
			a = 0;

			for (j = 0; j < n1; j++) {
				c = cos[a];
				s = sin[a];
				a += 1 << (m - i - 1);

				for (k = j; k < n; k = k + n2) {
					t1 = c * x[k + n1] - s * y[k + n1];
					t2 = s * x[k + n1] + c * y[k + n1];
					x[k + n1] = x[k] - t1;
					y[k + n1] = y[k] - t2;
					x[k] = x[k] + t1;
					y[k] = y[k] + t2;
				}
			}
		}
		Complex[] result = new Complex[n];
		for(i = 0; i < n; i++) {
			result[i] = new Complex(x[i] / Math.sqrt(n), y[i] / Math.sqrt(n));
		}
		return result;
	}
	
	/**
	 * 使用快速Fourier变换算法对复数向量作离散Fourier逆变换
	 * @param source 复数向量, 要求长度n必须为2的幂次, 即n = 2^k, k为整数
	 * @return 离散Fourier逆变换后的向量
	 */
	public static Complex[] ifft(Complex[] source) {
		int i, j, k, n1, n2, a;
		double c, s, t1, t2;
		int n = source.length;
		Pair<double[], double[]> twoArrays = Complex.toTwoArrays(source);
		double[] x = twoArrays.getFirst();
		double[] y = twoArrays.getSecond();
		
		int m = (int) (Math.log(n) / Math.log(2));

		// Make sure n is a power of 2
		if (n != (1 << m))
			throw new RuntimeException("FFT length must be power of 2");

		// precompute tables
		double[] cos = new double[n / 2];
		double[] sin = new double[n / 2];

		for (i = 0; i < n / 2; i++) {
			cos[i] = Math.cos(-2 * Math.PI * i / n);
			sin[i] = Math.sin(-2 * Math.PI * i / n);
		}
		// Bit-reverse
		j = 0;
		n2 = n / 2;
		for (i = 1; i < n - 1; i++) {
			n1 = n2;
			while (j >= n1) {
				j = j - n1;
				n1 = n1 / 2;
			}
			j = j + n1;

			if (i < j) {
				t1 = x[i];
				x[i] = x[j];
				x[j] = t1;
				t1 = y[i];
				y[i] = y[j];
				y[j] = t1;
			}
		}
		// FFT
		n1 = 0;
		n2 = 1;
		for (i = 0; i < m; i++) {
			n1 = n2;
			n2 = n2 + n2;
			a = 0;

			for (j = 0; j < n1; j++) {
				c = cos[a];
				s = sin[a];
				a += 1 << (m - i - 1);

				for (k = j; k < n; k = k + n2) {
					t1 = c * x[k + n1] + s * y[k + n1];
					t2 = -s * x[k + n1] + c * y[k + n1];
					x[k + n1] = x[k] - t1;
					y[k + n1] = y[k] - t2;
					x[k] = x[k] + t1;
					y[k] = y[k] + t2;
				}
			}
		}
		Complex[] result = new Complex[n];
		for(i = 0; i < n; i++) {
			result[i] = new Complex(x[i]/Math.sqrt(n), y[i]/Math.sqrt(n));
		}
		return result;
	}
	
	public static Complex[][] fft2(Complex[][] source) {
		int m = source.length;
		int n = source[0].length;
		Complex[][] result = new Complex[m][];
		for(int i = 0; i < m; i++) {
			result[i] = fft(source[i]);
		}
		for(int i = 0; i < n; i++) {
			Complex[] col = Complex.Col(result, i);
			Complex[] ffCol = fft(col);
			Complex.setColumn(result, ffCol, i);
		}
		return result;
	}

	public static Complex[][] ifft2(Complex[][] source) {
		int m = source.length;
		int n = source[0].length;
		Complex[][] result = new Complex[m][];
		for(int i = 0; i < m; i++) {
			result[i] = ifft(source[i]);
		}
		for(int i = 0; i < n; i++) {
			Complex[] col = Complex.Col(result, i);
			Complex[] ffCol = ifft(col);
			Complex.setColumn(result, ffCol, i);
		}
		for(int i = 0; i < m; i++) {
			for(int j = 0; j < n; j++) {
				result[i][j].real /= m * n;
				result[i][j].image /= m * n;
			}
		}
		return result;
	}
}
